Degrees & their distribution.

Degree of a vertex, in graph theory, refers to the number of edges that has the vertex as one of its ends. An in-degree and out-degree qualifies it, by counting only edges that start and end at the vertex. Degree distributions are important for graphs as they tell us about how connected the graph is and which if any of the vertices are more important than others.

degree(trips_graph, mode = 'out') %>%
  as.tibble %>%
  ggplot()+
  geom_density(aes(x=value))

degree.distribution(trips_graph) %>% head() # Check out what the degree.distribution produces and how to interpret the results.
# [1] 0.000000000 0.003880983 0.005174644 0.002587322 0.001293661 0.001293661

In particular, you want to look for the outliers in the distribution. For example to find the stations with only end trips but no start trips, i.e stations that are solely destinations

V(trips_graph)$name[degree(trips_graph, mode="out") == 0 & degree(trips_graph, mode="in") > 0]
#  [1] "3245" "3267" "3275" "3192" "3276" "3184" "3214" "3651" "3198" "3199"
# [11] "3203"

We can also visualise these stations that are popular origins on the map. However, we need to attach the values to the station locations tibble.

tmp <- degree(trips_graph, mode = 'out') %>%
  as.tibble() %>%
  rename(Outdegree = value)%>%
  mutate(stid = V(trips_graph)$name %>% as.numeric()) 


station_locs <- station_locs %>% 
  left_join(tmp, by='stid')

library(leaflet)

 Npal <- colorNumeric(
   palette = "Reds", n = 5,
   domain = station_locs$Outdegree
 )
 
station_locs %>%
    leaflet()  %>%
  addProviderTiles(providers$Stamen.TonerLite, group = "Basemap") %>%
   addCircles(
     lng = station_locs$lon,
     lat = station_locs$lat,
     radius = (station_locs$Outdegree - mean(station_locs$Outdegree))/sd(station_locs$Outdegree) * 30,
     fillOpacity = .6,
    fillColor = Npal(station_locs$Outdegree),
     group = 'Stations',
     stroke=FALSE
   ) %>%
  addLegend("topleft", pal = Npal, values = ~Outdegree,
            labFormat = function(type, cuts, p) {
              n = length(cuts) 
              paste0(prettyNum(cuts[-n], digits=0, big.mark = ",", scientific=F), " - ", prettyNum(cuts[-1], digits=0, big.mark=",", scientific=F))
            },
            title = "Out degree/Trip Starts",
            opacity = 1
  )

# Careful with division. If you have a nodes that are only origins, i.e. in-degree is 0, you will have a divison by 0 problem.

V(trips_graph)$name[degree(trips_graph, mode="out") / degree(trips_graph, mode="in") > 1.2]
#  [1] "216"  "258"  "266"  "289"  "339"  "356"  "391"  "397"  "399"  "443" 
# [11] "469"  "3152" "3161" "3164" "3177" "3302" "3316" "3366" "3536" "3539"
# [21] "3542" "3620" "3623"

Exercise

• Identify the popular origin stations of the trips by different times of the day and day of the week
• Modify the above code slightly to identify the popular destinations by different times of the day and day of the week.

Other centrality measures

station_locs$eigencentrality <- eigen_centrality(trips_graph)$vector

library(car)
k <- page_rank(trips_graph)$vector %>% Boxplot()

station_locs[k, 'stname']
# # A tibble: 10 x 1
#    stname                
#    <chr>                 
#  1 Pershing Square North 
#  2 West St & Chambers St 
#  3 E 17 St & Broadway    
#  4 Broadway & E 22 St    
#  5 W 21 St & 6 Ave       
#  6 12 Ave & W 40 St      
#  7 Kent Ave & N 7 St     
#  8 Broadway & E 14 St    
#  9 W 33 St & 7 Ave       
# 10 Central Park S & 6 Ave

k <- closeness(trips_graph, mode='all') %>% Boxplot()

station_locs[k, 'stname']
# # A tibble: 7 x 1
#   stname                
#   <chr>                 
# 1 Liberty Light Rail    
# 2 Newport Pkwy          
# 3 Essex Light Rail      
# 4 NYCBS DEPOT - DELANCEY
# 5 Morris Canal          
# 6 Columbus Drive        
# 7 Bike Mechanics at Riis

library(CINNA)
proper_centralities(trips_graph)
#  [1] "Alpha Centrality"                                
#  [2] "Bonacich power centralities of positions"        
#  [3] "Page Rank"                                       
#  [4] "Average Distance"                                
#  [5] "Barycenter Centrality"                           
#  [6] "BottleNeck Centrality"                           
#  [7] "Centroid value"                                  
#  [8] "Closeness Centrality (Freeman)"                  
#  [9] "ClusterRank"                                     
# [10] "Decay Centrality"                                
# [11] "Degree Centrality"                               
# [12] "Diffusion Degree"                                
# [13] "DMNC - Density of Maximum Neighborhood Component"
# [14] "Eccentricity Centrality"                         
# [15] "Harary Centrality"                               
# [16] "eigenvector centralities"                        
# [17] "K-core Decomposition"                            
# [18] "Geodesic K-Path Centrality"                      
# [19] "Katz Centrality (Katz Status Index)"             
# [20] "Kleinberg's authority centrality scores"         
# [21] "Kleinberg's hub centrality scores"               
# [22] "clustering coefficient"                          
# [23] "Lin Centrality"                                  
# [24] "Lobby Index (Centrality)"                        
# [25] "Markov Centrality"                               
# [26] "Radiality Centrality"                            
# [27] "Shortest-Paths Betweenness Centrality"           
# [28] "Current-Flow Closeness Centrality"               
# [29] "Closeness centrality (Latora)"                   
# [30] "Communicability Betweenness Centrality"          
# [31] "Community Centrality"                            
# [32] "Cross-Clique Connectivity"                       
# [33] "Entropy Centrality"                              
# [34] "EPC - Edge Percolated Component"                 
# [35] "Laplacian Centrality"                            
# [36] "Leverage Centrality"                             
# [37] "MNC - Maximum Neighborhood Component"            
# [38] "Hubbell Index"                                   
# [39] "Semi Local Centrality"                           
# [40] "Closeness Vitality"                              
# [41] "Residual Closeness Centrality"                   
# [42] "Stress Centrality"                               
# [43] "Load Centrality"                                 
# [44] "Flow Betweenness Centrality"                     
# [45] "Information Centrality"                          
# [46] "Dangalchev Closeness Centrality"                 
# [47] "Group Centrality"                                
# [48] "Harmonic Centrality"                             
# [49] "Local Bridging Centrality"                       
# [50] "Wiener Index Centrality"

Community Detection

Finding a community structure in a network is to identify clusters of nodes that are more strongly connected within the cluster than to nodes outside the cluster. There are many algorithms that identify community structures. Some work on weigthed networks and some work on digraphs. In any case it is better to convert multigraph (graph with two nodes having multiple edges) into a simple graph with edge weights, representing the number of trips.

E(trips_graph)$weight <- 1
station_graph  <- simplify(trips_graph, edge.attr.comb="sum")

ecount(station_graph) == ecount(trips_graph)
# [1] FALSE
sum(E(station_graph)$weight) == ecount(trips_graph)
# [1] TRUE
all.equal(vcount(station_graph), vcount(trips_graph))
# [1] TRUE

is.directed(station_graph)
# [1] TRUE

Walktrap algorithm finds the communities by random walks in the graph. The intuition is that random walks are more likely to be within a community that across communities.

wlktrp_mmbr <- data.frame (clstrmm  = cluster_walktrap(station_graph)$membership %>% as.factor, 
                          stid = V(station_graph)$name %>% as.numeric()
                          ) %>% as.tibble()


wlktrp_mmbr$clstrmm %>% summary()
#   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15 
#  80 218 207  61   3 195   1   1   1   1   1   1   1   1   1

It seems like there are 5 distinct and non-trivial clusters of stations. You can visualise these results by joining them back to tibbles with location information.

station_locs <- station_locs %>%
                inner_join(wlktrp_mmbr, by='stid')

library(ggmap)
nybg <- get_map(location=c(-74.0060,40.7128), maptype = 'toner-lite', zoom = 11)

library(RColorBrewer)
numColors <- levels(station_locs$clstrmm) %>% length()
myColors <- colorRampPalette(brewer.pal(8,"Dark2"))(numColors) # Expanding the number of colors available from 8 

  ggmap(nybg)+
  geom_point(aes(x=lon, y=lat, color=clstrmm), alpha=.9, data=station_locs)+
  scale_color_manual(values = myColors)+ 
     scale_x_continuous("", breaks=NULL)+
   scale_y_continuous("", breaks=NULL)+
   theme(panel.background = element_rect(fill='white',colour='white'), legend.position = "none") +
  labs(colour = "Communities") 

Exercise

• Do the clusters change if you use a different algorithm?
• Are the clusters different for customers and subscribers? Men and Women?

Unemployment rate

Unemployment rate is trickier, because it needs to track both the labour force as well as persons unemployed. It turns out that ACS only provides this data for males and females separtely and by differnt age categories. So we need to assemble them and combine them.

# variables for labour force for age groups 16-64
lf_m <- paste("B23001_", formatC(seq(4,67,7), width=3, flag="0"), "E", sep="")
lf_f <- paste("B23001_", formatC(seq(90,153,7), width=3, flag="0"), "E", sep="")

lf <- reduce(
  map(NC_SC, function(x){
    get_acs(geography='tract', variables = c(lf_m, lf_f), state=x, year=2016)
  }),
  rbind
)

# variables for unemployed for age groups 16-64
unemp_m <- paste("B23001_", formatC(seq(8,71,7), width=3, flag="0"), "E", sep="")
unemp_f <- paste("B23001_", formatC(seq(94,157,7), width=3, flag="0"), "E", sep="")

unemp <- reduce(
  map(NC_SC, function(x){
    get_acs(geography='tract', variables = c(unemp_m, unemp_f), state=x, year=2016)
  }),
  rbind
)

lf_t <- lf %>% 
  group_by(GEOID) %>%
  summarize(lf_est = sum(estimate, na.rm=T))

unemp_t <- unemp %>% 
  group_by(GEOID) %>%
  summarize(unemp_est = sum(estimate, na.rm=T))

unemp_rate <- merge(lf_t, unemp_t, by='GEOID') %>% 
  filter(lf_est >0) %>%
  mutate(unemp_rate = unemp_est/lf_est)

Exercise

• Compute unemployment rate for more meaningful classes (Women, Youth etc.) and visualise it. Are there differences in spatial distribution?

We can arbitarily define distress to be a combination of poverty and unemployment. In this instance, I am picking 20% and 10% as the threshold respectively.

tract_stats <- merge(pov_rate, unemp_rate, by='GEOID')
distressed_tracts <- tract_stats %>%
                     filter(unemp_rate > .10 & pov_rate > .20)

mapview(distressed_tracts)

There are reasons to presume that clusters of distress is more problematic than isolated distress. We could use the contiguity of distressed tracts to our advantage to figure where these clusters are and how they are shaped.

It turns out that the contiguity calculations in sf are still a little slower than older packages such as spdep. So I will convert the simple features into a Spatial object and use older packages.

library(rgeos)
library(spdep)

dt <- as(distressed_tracts, "Spatial") # Convert to spatial
temp1 <- gUnarySTRtreeQuery(dt) # Construct list of polygons that are likely to intersect/touch by first looking at bounding boxes.
dt_nb <- poly2nb(dt, queen=FALSE, foundInBox=temp1, row.names = dt$GEOID) # Construct a neighborhood object

plot(dt)

plot(dt_nb, coordinates(dt), col='red', points=F)

distressed_tracts_graph <- dt_nb %>%  
  nb2mat(zero.policy=TRUE, style="B") %>%
  graph.adjacency(mode='undirected', add.rownames=NULL) # construct an undirected graph from the adjacency matrix.

cl <- clusters(distressed_tracts_graph)

ggplot() + 
  geom_histogram(aes(cl$csize))

sum(dt$GEOID!= names(cl$membership)) # this is 0. So we can safely assign the clustermembersip variable by rows rather than doing a merge.
# [1] 0
dt$clustermembership <- cl$membership

Clusters function decomposes the graph in subgraphs, i.e. nodes belong to a subgraph, if there is a path between them. If there is not, they belong to different subgraphs.

pal <- colorRampPalette(brewer.pal(8, "Dark2"))
numColors <- cl$membership %>% unique() %>% length()

dt %>% st_as_sf() %>%
  mapview(zcol="clustermembership", 
          col.regions=pal(numColors),
          legend = FALSE,
          lwd = 0
          )

Wihtin in each component of the graph, you can identify a community structure if you wish, using techniques from above.

In some instances, it may be useful to see if the clusters are stringy’ or if they are ‘blobs’. i.e., if the clusters follow linear features (such as roads, valleys, rivers etc) or if they adjacency is shared among multiple polygons of the same component. To do this, we could rely on diameter of the graph, which is the longest path between any two nodes in graph. If the graph is stringy, it will have a large diameter.

In this instance, we want to compute the diameters of each component separately. (diameter of the graph is \(\infty\), as it is not fully connected. )

diameters_of_graphs <-   distressed_tracts_graph %>%
       decompose() %>%
       sapply(function(x)diameter(x)) %>% 
       as.tibble()

names(diameters_of_graphs)[1] <- 'dia'

diameters_of_graphs$cl_membership <- 1:(distressed_tracts_graph %>% decompose() %>% length())

dt <- merge(dt, diameters_of_graphs, by.x='clustermembership', by.y='cl_membership')

plot(dt[dt$dia > 10,'GEOID'], col='red') # Quickly showing the stringy and blobby (real words) distress
plot(dt[dt$dia < 3,'GEOID'], col = 'blue', add = T)

You will notice that most of large spatial concentrations of distress is primarily in the rural ares (midlands of South Carolina and Eastern Carolina Plains). However, this may be misleading. There are concentrations of distress within the urban regions. Since census tracts are very small, they are more likely to be clustered and separated. So the unequal sizes of geography affects the conclusions you can draw.

Exercise

• Limiting your analysis to Metropolitan Areas in Carolinas, what conclusions can you draw about the shape of the clusters?

• Does limiting your analysis to a single state, affect your results?

Instead of categorising the areas of distress, one could use local Moran's I to identify locations high values of a continous variable (say poverty rate) surrounded by other high values areas (or low value areas) and test their statistical signifance. local.moran is a function in spdep and relies on the neighbourhood graph. I leave this as an exercise.